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Abstract 

Determining if two histograms are consistent, whether they have been drawn from 
the same underlying distribution or not, is a common problem in physics. Existing 
approaches are not only limited in power but also inapplicable to histograms filled 
with importance weights, a common feature of Monte Carlo simulations. From a 
Bayesian perspective, the comparison between a single underlying distribution and 
two underlying distributions is readily solved within the context of model comparison. 
I introduce an implementation of Bayesian model comparison to the problem, including 
the extension to importance sampling. 

1 Histogram Comparison 

The histogram is a fundamental tool in physics, providing a powerful yet accessible means of 
non-parameteric inference. Often, however, analyses are concerned not with the underlying 
distribution itself but rather the consistency between multiple histograms. Comparing dif- 
ferent analysis techniques, for example, requires the comparison of each respective analysis 
output given the same input data. Perhaps the most ubiquitous example is the validation 
of Monte Carlo simulations against data. 

Although numerous approaches to the problem exist, none are without their limitations. 

2 Orthodox Approaches 

2.1 The Komolgorov-Smirnov Test 

When comparing the consistency of two one-dimensional samples, yi and y~2, standard fre- 
quentist methodologies turn to significance testing, in particular the Kolmogorov-Smirnov 
(KS) test [I]. Here the null hypothesis that the two samples are consistent with draws 
from the same distribution is accepted or rejected based on the sampling distribution of 
the test statistic 

q = max |F yi (x) - F y2 (x) \ , 
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where F z (x) is the empirical distribution of the sample z, 



n 




with I (xi) the indicator function 




0, x < Xi 

1 •j CC '^i, 



In practical applications, the sampling distribution for q is usually approximated with 
a large n limit such that the KS test becomes extremely conservative (in more technical 
terms, the Type II error of the test is small only when the samples are large). Moreover, 
the sampling distribution of the weighted samples from importance sampling cannot be 
approximated so easily and the proper application of the KS test to simulated events is 
infeasible. 

Given histograms instead of the individual samples, the empirical distributions must 
be approximated by the bin contents themselves and the test statistic q reduces to the 
maximum difference between bins. Information in all other bins is ignored, no matter 
now relevant it might be. Considering also the usual faults of frequentist significance tests 
IU Hj, the application of the KS test to the problem of histogram comparison leaves 
much to be desired. 

2.2 Other Approaches 

Various procedures have also been developed in various communities. One approach sub- 
tracts one histogram from the other before performing a significance test on the hypothesis 
of a constant residual across all bins. The sampling distribution of the residuals, however, 
does not admit familiar tests such as x 2 outside of the limit of large bin contents. Still, 
the limit is often assumed and the resulting x 2 taken as a measure of consistency between 
the two histograms. 

A similar approach utilizing the quotient of the two histograms fares even worse. The 
usual linearized Gaussian approximations to the resulting bin quotients are not very accu- 
rate, and the subsequent x 2 test is even more misleading than for the bin residuals. More 
accurate considerations [5] become considerably more involved, and significance testing 
much more challenging. 

Both approaches are further limited by the dependence on frequentist significance test- 
ing. 
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3 The Bayesian Perspective 



From the Bayesian perspective, the problem of histogram comparison becomes one of model 
comparison between 5, the model where the two histograms are drawn from the same 
distribution, and S, the model where they are drawn from distinct distributions. The 
probability of S follows from an application of Bayes' Theorem on the model evidences, 
p (m, n\S) and p (n, m|<S) , 

/ C i \ p(m,n\S)p(S) 
p (o|m, n) - 



p (m, n\S) p{S)+p (m, n\S) p (S) ' 

where m and n are the bin populations of the two histograms. In practice, a uniform prior 
is taken and the probability reduces to 

/ci ^ p(m,n\S) 
p(«S|m,n) ~ 



p (m, n|<S) + p (m, n|<S) 

Subsequently, S best models the data when p(S\m, n) > 0.5 while S is superior when 
p(S\m, n) < 0.5. When p(S\m, n) = 0.5, however, the interpretation is more subtle. The 
intermediate case can arise not only when each bin is independently ignorant, but also 
when individual bins are best described by different models and the discrepancies negate 
the respective contributions to the model evidences. 

In order to discriminate between ignorance and disagreement, first note that the un- 
derlying multinomial distributions modeling the histogram bin contents are well described 
by independent Poisson distributions. The resulting evidences conveniently factor, 

N 

p(m,n\M) = Y[p(mi,ni\M) , 

i=i 

which enables an explicit mixture model, 



N 

p(m,n|vr) = JJp (m*, 7^71-) 
i=i 
N 

p(m,n|7r) = [up (m^n^S) + (1 - n)p (m^n^S)] , 



i=i 



with the mixture posterior 

p (7r|m, n) oc p (m, n|7r) p (it) . 



3 




(a) (b) (c) 

Figure 2: When all bins agree the posterior is either (a) concentrated towards tt = and 
favoring S , (b) indifferent, or (c) concentrated towards tt = 1 and favoring S . 



In the case of a single bin and a uniform prior, the mixture posterior is linear in tt and 
the possible posteriors fall into three classes. Those with a negative slope are best modeled 
by S (Fig[Th.), those uniform across tt are indifferent between the two models (Fig [lb), and 
those with a positive slope favor S (Fig Ht) . 

In the general case, each bin is modeled independently. If all bins favor S then the 
combined product quickly concentrates towards small tt with a mode exactly at tt = (Fig 
[2^,). Conversely, full support of S produces a posterior concentrated towards tt = 1 with a 
mode exactly at the boundary (Fig [2b). Any bins indifferent between the models contribute 
little to the final shape, unless all the bins are ignorant in which case the posterior remains 
uniform (Fig [2b). When bins are best modeled by different models, however, the posterior 
becomes peaked away from the boundaries with the exact location of the mode indicating 
the overall superior model (Fig [3]) . 

Diagnosing the best model is then straightforward. When a mode arises at the bound- 
aries, the bins are consistently described by a single model: S for tt = and S for tt = 1. 
A mode away from the boundaries indicates mixed agreement, implying that the two his- 
tograms are not entirely consistent. 
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Figure 3: A local mode arises when individual bins are best described by different models. 



Note that the non-mixture posterior can be recovered by comparing the boundaries of 
the mixture posterior, 

/ c , \ p(7r = l|m,n) 

p (o|m, nj 



p(ir = 0|m, n) + p (ir = l|m, n) 

The additional information contained in the shape of mixture posterior underlies the 
benefit of the extension to a mixture model. Not only does the mixture model break the 
the degeneracy between ignorance and balancing disagreement, but it can also reveal more 
subtle disagreements that might otherwise be obscured in the model posterior. 

3.1 Comparing Data 

Actually computing the mixture posterior requires the bin by bin evidences, p(mi,rii\S) 
and p (mi,rii\S) . 

3.1.1 The Same Source Model Evidence 

When both bin populations are drawn from the same distribution, the joint likelihood is 
given by 



p(rrii,ni\\,S) = p (mi\\)p (rii\X) 

A m, e -A \n l(} -\ 



m»! nil 

yrii+rii e -2\ 

mj.riil 
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Assigning a uniform prior, 

_ / (/3- a) -1 , a < X < p 



p(X\S) 
yields the posterior 



p(\\m.i,ni,S) oc p(mi,rii\X,S)p(\\S) 

13 A \ \m,+n, p -2X 



0, else 



fdAA 

J a 



2 -(m i+ n;+i) r ( m . +ni + l) [7 ( m - + n - + l, 2/3) - 7 (mi + n { + 1, 2a)] ! 
where 7 is the normalized lower incomplete gamma function, 



7 j«, dttn _ le _ t - 
The evidence follows, 



/ , C N p(m i ,nj|A,<S)p(A,5) 
p(mj,nj|<S) = 



p(X\mi,ni,S) 

1 2-( m » +n ' +1 ) 7 (mj + «i + 1, 2/3) - 7 (nii + rc, + 1, 2a) 
(3-ami + rii + l Be (m^ + 1, nj + 1) 



where Be (x, y) is the Beta function, 

Be{x,y) = C dtt x {\ -tf- y . 
Jo 

In practice the prior support, a < X < (3, is chosen to encompass the bulk of the likeli- 
hood. Making the support as small as possible while keeping the posterior normalization 
practically constant implements a crude form of model selection that improves the model 
without limiting its power. 



6 



3.1.2 The Two Source Model Evidence 

Given two distinct sources, the joint likelihood factors into 



p(mi,rii\\,S) =p(mi\X) -p(ni\p) 
rriil nil 

The two sources require two independent priors, 

p{X,p\S) =p(\\S) -p{p\S) 

(S - 7)" 1 ■ (C - e)' 1 , J<\<5,e<u<( 
0, else 

and the joint posterior becomes 

p (X, u\nn,ni,S) = p (\\mi,m,S) ■ p (u\mi,rii,S) 

ccp (m i ,n i \X,S)p (X\S) ■ p (m,, n^p, S) p (p\S) 

d\ X mi e~ x dp, n n ie-» 

A m le -A p n i e ~^ 



T ( mi + 1) [7 (rrn + 1, 5) - 7 (m< + 1, 7)] r (n* + 1) [7 K + 1, C) - 7 + 1, 0] ' 
The evidence is then 



/ _ p{m h ni\X,pS)p(X,pS) 



p(X,u\mi,ni,S) 
p(mi,m\X,S)p (AS) p{m il n i \u,S)p (pS) 

p(X\nii,ni,S) p(u\mi,m,S) 
T (mj + 1) [7 (im + 1, 6) - 7 (mj + 1, 7)] T (m + 1) [7 (n< + 1, Q - 7 (»» + 1, e) 

r ( mi + 1) (<5 - 7) r (ni + 1) (C - e) 

[7 (m^ + 1, 8) - 7 (mj + 1, 7)] [7 (rij + 1, C) - 7 K + 1, <Q] 
(«5- 7 ) (C-c) 
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3.2 Importance Sampling 

3.2.1 Sampling as Monte Carlo Estimation 

Cross sections defining the physics from which the data are generated are equivalent to the 
probability distribution 

. 1 da (x) 

P(x) = -j , 

a ax 

where x defines the phase space of the given physics. 

Given an integrated luminosity C, the expected number of events falling into a particular 
histogram bin is then 



A= / dxC da{x) 



-L 



dx 

f b ~ 
= dx Cap (x) 

J a 

dxl ab (x) Cap(x) 

> 

dxN (x)p (x) , 

J — oo 

where 



'OO 

oo 



lab 0) 



1, a < x < b 
0, else 



Now consider a Monte Carlo simulation generating n = Ca events from the differential 
cross section da/dx. The samples admit the Monte Carlo estimate of the expected counts, 



N = l±N(x t ) 

i=l 
1 H 

n ^— ' 

i=i 

n 

= Iaf > ( Xi "> ' 



i=i 

which is just the number of events falling into the bin a < x < b. 
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The estimate is unbiased, 



E p [N] = A, 



with the variance 



Var p [iV] = -Var p [N (x)\ 



= -(e p [N 2 (x)}-E p [N(x)] 2 ) 
= l(fj 2 ab (x) (Ca) 2 p(x)-X 2 




Ca Iab(x) Cap(x) - A 2 



— oo 



= I (nA - A 2 ) 

n 

= A - A 2 /n. 



In the large n limit, the second term becomes negligible and the sampling distribution 
of the Monte Carlo estimate converges to a Poisson distribution, 



Note that this limit does not require large bin contents. In fact, the large n limit is exactly 
the same rare-events limit in which the multinomial distribution reduces to independent 
Poisson distributions. The Monte Carlo samples, then, are entirely analogous to the data 
and inferring the underlying rate A is the same in both cases. 

Importance sampling is slightly more subtle. Here events are not drawn from the 
physics distribution p (x) but instead an auxiliary distribution q (x). The expected number 
of events becomes 
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with the resulting Monte Carlo estimate 



N = -i2N(x l )w(x i ) 
t=\ 

1 n 

= - ^2CaI ab (xi)w(xi) 
n i=i 

La n 

= {Xi)w{Xi) 

i=l 

n 

= ^2 Iab ( Xi ~) w ^ ■ 
i=i 

Note that this estimate is just the sum of the weights of events falling into the chosen bin. 
As above, the importance sampling estimate is an unbiased estimator, 

E q [N] = A, 

but with the variance 



VarJiV] = ^Var, [N (x) w (x) 



X - (E, [N 2 (x) w 2 (x)] - E q [N (x) w (x)} 2 ) 

i ^La J I ab (x) La w (x) p (x) - A 2 ^ 
- J I ab (x) La w (x) p (x) - A 2 ^ 



n 



oo 

I a b (x)Law (x)p(x) — X 2 /n. 



If the weighting function w (x) is relatively constant across the bin a < x < b then 
w (x) can be approximated by its ensemble average, 

^ _ Yli=l tab W (Xi) 



Y2i=l lab ( x i 
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and the remaining integral becomes 



/OO POO 
I ab (x) Ca w (x) p (x) ^ J I ah {x)Cawp{x) 
-oo J — oo 

/oo 
I ab (x) Cap(x) 
-oo 

ps u>A. 

The variance is then 

Var g [iV] RiwA - A 2 /n, 
and in the large n limit reduces to 

Var 9 [iV] pa u)A. 

Note that if the weights are unity, in which case the importance sampling reduces to stan- 
dard Monte Carlo sampling, then the variance reduces to the previous result as consistency 
would demand. 

The resulting sampling distribution is approximately Gaussian, 

p(n\X) = . 1 exp 

The above sampling distribution admits a Bayesian approach to simulated histograms 
and consequently the comparison of simulated histograms. Computation of the evidence 
ratio follows as in the case of data comparison. 

3.3 Comparing Data to Simulation 

The importance sampling distribution for the simulated histogram bins requires not only 
the sum of the weights falling into that bin, 

n 

rii = ^2l ai bi {Xj)w(xj), 

3=1 

but also the ensemble average of the weights, 

_ _ TTj^hjbi (xj)w(xj) 

2^j=l 1 a i b i { x j) 

As before, the sampling distribution of the data histograms bins rrii is taken to be 
Poisson. 



l(iV-A) 
2 wX 
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3.3.1 The Same Source Model Evidence 

Here the joint likelihood becomes 



p(mi,rn\X,S) = p (mi\\) p (ni\\) 



nu „— A 



A m »e 



A uniform prior, 



rn.il y/2irwiX 



(/3-a)" 1 , a<X</3 
0, else 



exp 



1 (ni - xy 

2 WiX 



P (X\S) 



gives the posterior 



p(\\m,i,ni,S) oc p(mi,rii\X,S)p(X\S) 





A - 


1 (rij— A) 2 " 

2 WiX 




dAA m '-5exp 


\ l(«.-A) r 
A 2 ?D; A 



resulting in the evidence 
p(mi,rii\S) 



1 



(/3 — a) \/2nwi 



dX X mi 2 exp 



A- 



1 (nj - A)- 

2 u),;A 



Note that limiting the prior support not only improves the model but also eases the 
computational demands of the required numerical integration. 

3.3.2 The Two Source Model Evidence 

Given two distinct sources, the joint likelihood factors into 



p{m h rii\X,S) =p(rrii\X) -p(rii\n) 
A m 'e" A 1 

mil y/2-KWiX 



exp 



1 (ni - xy 

2 WiX 
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The two sources require two independent priors, 



p(A,/*S) =p(X\S) .p(»\S) 

f (S-^-iC-ey 1 , 7 <A<<5,e</x<C 
\ 0, else 

and the joint posterior becomes 



p(X,fjb\mi,rii,S) = p(\\nii,ni,S) ■ p (fi\mi,rii,S) 

oc p (mi,m\X,S) p (X\S) -p (mi,rii\fx,S) p (/z|<S) 



1 

\mi e -\ jJL 2 exp 


1 {rii-jj,) 2 

2 uiifi 




f 7 5 dAA m *e" A j^d^-iexp 


1 (rii-n) 2 

2 Wifj, 



l 

yni e -\ /J, 2 exp 


1 (rii-n) 2 

2 Win 




T(mi + 1) [7 (mi + 1,5) - 7 (m ; + 1,7)] J^^^-h exp 


1 (ni-fi) 2 

2 Wifi 



yielding the evidence 



t icA _v{™i,n%\X,\i,S)p{\,\i,S') 



p (X,fi\mi,m,S) 
p (mi,rii\X,S)p(X\S) p (mi,rii\fi,S) p (fi\S) 
p(X\mi,rii,S) p(fi\mi,rii,S) 



T (mi + 1) [7 (mi + 1,5) - 7 (mi + 1,7)] ■ y d/x/x 2 exp 



3.4 Comparing Simulation to Simulation 

Comparing two simulated histogram bins requires the sum of the weights and ensemble 
average for both the first bin, 

n 

mi = ^I a . b . ( Xj )w(xj) 
3=1 

TI l j=l I a i b i (Xj)w(Xj) 



13 



and the second bin, 



ni = ^2la t b t (Xj)v(Xj) 
3=1 



3.4.1 The Same Source Model Evidence 

In this case the joint likelihood becomes 



p(mi,rii\X,S) = p(mi\\)p(rn\\) 



1 



A uniform prior, 



1 



p(X\S) 



cxp 



1 (rrii - \y 

2 WjX 



\j2irviX 



cxp 



1 (n t - xy 

2 ViX 



cxp 



1 ( (mi - X) (rii - X)' 



2A \ Wi 



v i 



((3-a)~ l , a<X<f3 
0, else 



gives the posterior 



p(X\mi,rii,S) <xp(mi,rii\X,S)p(X\S) 



X 1 exp 




( '(mi-X) 2 (rn-X) 2 \ 

^ Wi Vi J 






£fdAA-iexp 




1 /(m t -A) 2 K-A) 2 \l 
2A y Wi Vi J 



with the evidence 



1 r - 

p(mi,rii\S) = — — j== I dAA 1 exp 

2vr (J3 - a) y/mvi J a 



J_ / (mi - A) 2 (ni - X) 

2 A I iDi 
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3.4.2 The Two Source Model Evidence 

With two distinct sources the joint likelihood factors into 



p(mi,m\X,S) =p(m;|A) -p(ni\p) 



1 



\ / 2irWi\ 
1 



exp 



1 (nii - \y 



1 



2-Ky/WiV. 

Taking two independent priors, 



A 2 exp 



2 WiX 

1 (m, - A) 



\J2-KVill 
2 



exp 



1 (rij - /j,y 

2 Uj// 



2 If; A 



1 


1 (n, - n) 2 


• jU 2 exp 




2 Vi/x 



p(A,/,|5)=p(A|5)-p(/i|5) 

_/ (S-^-iC-e)' 1 , ^<X<6,e<p<( 

0, else 



gives the posterior 



p(X,fj,\mi,ni,S) = p(X\mi,rii,S) ■ p (fi\mi, rii,S) 

oc p (rrtj, nj| A, <S) p (X\S) ■ p (m*, n^/i, S) p (p\S) 



X 2 exp 


l(mi-A) 2 " 


l 

/i 2 exp 


1 (n—^t) 2 

2 Sj/i 




J 7 5 dA A 2 exp 


l(mi-A) 2 " 
2 WiX 


du /i 2 exp 


1 (n—^t) 2 

2 Si/i 



along with the evidence 

p(m i ,n i \X,fi,S)p(X,fi,S) 



p(m h rii\S) 



p (X,p\mi,rii,S) 
p (mi,ni\X,S) p (X\S) p (mi,n;|u,S) p (p\S) 



p (X\mi,ni,S 
((5 _ 7 )l. (C _ e )-' r* 



/ dA A 2 exp 



p (fi\m,i,ni,S) 

1 (mi - A) 2 



2 w,A 



r 



d/i jU 2 exp 
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4 Examples 



| input Histograms | | p{n > 1/2) = 0-475, = 0.497, p(S| m, n) = 




Figure 4: Comparison of two histograms drawn from different distributions. 



| Input Histograms | | Pfr > 1/2) = 0.999, ;r Mflp = 1, p(S| m, n) = i~| 




7t 



Figure 5: Comparison of two histograms drawn from the same distribution, with the red 
sample generated from importance sampling. 



| Input Histograms) 
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p(ir>1/2) = 0,J Mftp = 0,p(S|m, n) = | 
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Figure 6: Comparison of two histograms drawn from the same distribution, with both 
samples generated from importance sampling. 
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5 Model Ensemble Tests 



The capability of the Bayesian approach is evident when considering an ensemble of his- 
togram comparisons models. 

At each iteration in the ensemble a true model is randomly selected between S and 
S. In the former case, the two histograms to be compared are drawn from a randomly 
generated multinomial distribution with 50 bins, 

Pi 

P; - 



Pi ~U(0,1). 

The latter case samples one histogram from the randomly generated multinomial dis- 
tribution, but draws the second from a perturbed distribution with 

Qi 



* = ( i+ H r< -^ 

ri~U(0,l). 

The two distributions are kept similar in order to assess the performance of each algorithm 
when comparison is highly nontrivial. 

The two histograms generated for each model are then evaluated with three tests: 

Non-Mixture Posterior : p(S\m, n) > 0.5 

7TMAP > 0.97 

A 95% significance KS test 
implemented in the ROOT [6 J physics library 

False acceptances and false rejections were tabulated for four ensembles, each with 2000 
independent models but different total bin contents n. The respective rate posteriors were 
calculated assuming a binomial likelihood and beta prior, Be (0.5, 0.5); the posterior modes 
and 68.3% confidence intervals for each algorithm are plotted in Figure [71 

Note that in the two cases with the smallest statistics, n = 5, 000 and n = 10, 000 the 
variance of the bin contents equals or surpasses the expected difference between the two 
underlying distributions and the large false accept rates are to be expected. 

In all cases the Bayesian mixture test, vtmap > 0.97, is superior to KS. The non-mixture 
test outperforms KS in almost all cases as well, failing only in the small statistics ensembles 
where small discrepancies tend to be washed out by the bin content variances. 



Mixture MAP 
Kolmogorov-Smirnov 
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Figure 7: (a) False reject and false accept rate for each comparison algorithm. 
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6 Conclusion 



A Bayesian mixture model has been developed to test whether two histograms are con- 
sistent, in that they are more likely to have been drawn from a single distribution rather 
than two distinct distributions. The model is extended to handle histograms generated 
from importance sampling, resulting in a robust and powerful approach to the comparison 
of histograms populated by both data and simulation. Said power is evident with studies 
of a large model ensemble. 

A C++ implementation of the Bayesian mixture model utilizing the ROOT [6j data 
analysis framework is available at |http : //web . mit . edu/ -be t an/www/ code . html, 

7 Acknowledgements 

I thank Chris Jones, Steve Voinea, and Matt Walker for helpful discussion and comments. 

References 

[1] Papoulis, A. and Pillai, S. U. (2002) Probability , Random Variables, and Stochastic 
Variables. McGraw-Hill, New York 

[2] Bernardo, J. M. and Smith, A. F. M. (2000) Bayesian Theory. John Wiley & Sons, 
Inc., New York 

[3] Jaynes, E. T. (2003) Probability Theory: The Logic of Science, Cambridge University 
Press, New York 

[4] MacKay, D. J. C. (2003) Information Theory, Inference, and Learning Algorithms. 
Cambridge University Press, New York 

[5] Hinkley, D. V. (1969) Biometrika 56(3) 635-639 

[6] Brun, R. and Radermakers, F. (1997) Nucl. Instr. & Meth. A 389(1-2) 81-86. See also 
htt p : //root . cern . ch/[ 



19 



